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Evolution of density and metric perturbations in the background of high frequency oscillations of 
curvature in F{R) gravity is considered. In addition to the usual Jeans-like instability new effects of 
amplification of perturbations, associated with parametric resonance and antifriction phenomena, 
are found. 


I. INTRODUCTION 


Gravity modifications at large distances have been suggested [1] for an explanation of the accelerated cosmological 
expansion observed in the present-day universe. The idea is to add a non-linear function of curvature scalar, F{R), 
to the standard Einstein-Hilbert action: 

4rav = - ^/ [R + Fm , (1.1) 

such that the modified Einstein equations have an accelerated De Sitter-like solution in the absence of matter^. These 
theories contain an additional massive scalar degree of freedom (dubbed “scalaron”) beside the usual massless, spin-2 
graviton. In the first works [1] the function F{R) was taken in the form F{R) = —pL^/R, where /r is a constant 
parameter with dimensions of a mass and of the same order of magnitude as the present Hubble parameter (or inverse 
universe age), i.e. ^ \Rc\ ^ ^/tu, where tu fv 14 Gyr is the universe age. It was found [2], however, that in 
astronomical systems with even slightly varying mass density a very strong instability would develop, resulting in 
explosive solutions evidently incompatible with observations. 

To cure this shortcoming further modifications of GR have been suggested [3-5], for a review see Ref. [6]. Let us 
take as a guiding example the specific F{R) of Ref. [3]: 


F{R) = -XRc 



( 1 . 2 ) 


where Rc is a constant parameter with dimension of curvature and close by magnitude to the present cosmological 
curvature, A is a dimensionless constant of order unity and the power n is usually taken to be an integer (though 
not necessarily so). This choice of F{R) leads to a further problem, namely the solution of the equation of motion 
for the gravitational field must be singular with i? —> oo in the past to produce a reasonable late-time cosmology [7]. 
Similarly, it has been found that systems with rising mass/energy density will evolve to a i? —>■ oo singularity in the 
future [8]. The singularity can be avoided if one adds an extra term F{R) — F{R) — R?/{Qm?) [7]. At small R 
the system tends to evolve to higher values of curvature, but as \R\ grows the i?^-term eventually becomes dominant 
and pushes the system back to lower curvatures. This results in oscillating solutions R{t), possibly with very large 
amplitude. In our papers [9] we have found such oscillating solutions for the particular example of F{R) given by 
Eq. (1.2) with the addition of the R^ term. The properties of different versions of F{R) are such that they induce 
effective potentials for R having minimum at R equal to its GR value, R = RgR: and rising in both directions of lower 
and higher R. So an oscillating behaviour of i? is a generic feature of these modified gravity theories. We assume that 
the form of such oscillations is arbitrary, keeping their amplitude and frequency as free parameters. 

Gravitational instability in modified gravity has been considered in Refs. [10]. In our paper [11] we studied the Jeans 
instability in classical and modified gravity in a background with rising energy density, assuming that the background 
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metric changes slowly as a function of space and time compared to the typical frequency and Compton wavelength 
of the scalaron field. Here we investigate gravitational instability in a quickly oscillating curvature background. The 
time evolution of first-order perturbations is governed by a fourth-order differential equation (instead of the usual 
second order one), therefore new types of unstable solutions are expected, and have in fact been found. 

There appear not only the usual Jeans solution with a slightly modified (reduced) length scale, but also parametric 
resonance amplification of density perturbations, and an amplification of the perturbations due to the ’’antifriction” 
behaviour of the coefficients of odd derivatives (first and third) in the equation. 

The paper is organized as follows. In Section II we give a general description of density perturbations in a quickly 
oscillating background; in Section III we discuss the Jeans-like instability in modified gravity. In Section IV parametric 
resonance amplification of fluctuations is studied for harmonic and spiky oscillations of the background curvature. In 
Section V we investigate the effects of the antifriction induced by the change in the signs of the coefficients in front 
of the odd derivatives in the equation of motion. Finally, Section VI is dedicated to the discussion of the results 
and possible implications. In Appendix A we present the expressions for the metric, Christoffel symbols, Ricci and 
Einstein tensors for a Schwarzschild-like isotropic and homogeneous solution in modified gravity. 


II. DENSITY PERTURBATIONS OVER A QUICKLY VARYING BACKGROUND 

A. Basic framework and equations 


The modified Einstein field equations, derived from the action (I.l), are: 


(1 + F'p.) 


{R + F) {g^vF>aF)°‘ — — 


SttT, 




= T„ 


'-PI 


( 2 . 1 ) 


where F'j^ = dF/dR, is the covariant derivative and is the energy-momentum tensor of matter. 

We assume that the background metric weakly deviates from the Minkowsky one, while derivatives of the metric 
may be far from their GR values. In particular, R may be very much different from Rgr = ~F. Please note that 
these two conditions are only apparently contradictory; in fact R ~ d^g so even if the metric is approximately the 
Minkowsky one, the curvature may be very large in comparison with Rgr, e.g if the frequency of the oscillations of 
g and R is high, as in the case under scrutiny. We consider astronomical systems with |i?c| ^ |R| 'C m?. Both limits 
are natural for relatively dense systems, that is denser than the average cosmological background but much less dense 
than rn?pirn?. Note that m > 10^ GeV, see e.g. [12]. In this limit we have 

|F(i?)|«|R|, |F'(i?)|«l. (2.2) 


This is surely fulfilled for the model (1.2), for which at i? ^ Re- 


F{R) 


—XRc 


1 - 



R^ 

drrF 


Therefore, equation (2.1) can be approximated as: 

Gfiu -f — gfj,uD^)R = , 

where = R^i, — g^^R/2 is the Einstein tensor and 


= -^F'^r- 


(2.3) 


(2.4) 


(2.5) 


Once written in this form, the equation is largely independent of the specific model considered except of course for 
the value of w, and provided that along the background solution w ~ const., or rather io/uP' <C 1, and that the 
conditions (2.2) are fulfilled. In the example (2.3): 


1 6An(2n-fl) 

\r) 


( 2 . 6 ) 


Gosmological perturbations started to rise at the onset of the matter dominated epoch at redshift Zgq — 10"^ when 
Req/Rc ^ 10^^. So for m = 10^ GeV, which is its lower limit [12], w may be treated as practically constant if n > 3. 
For systems with energy density of about 1 g/cvep the frequency uj would remain constant even for n > 1. If w rises 
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with time, the perturbations would rise even faster than obtained below. If (2.6) is dominated by the second term in 
square brackets, that is 

^6An(2n + l) [Jij ’ 

then the frequency might depend crucially on time and the approximation of constant frequency would be invalid, 
although |i?| ^ |i?GR| makes the second term in Eq. (2.6) smaller. This corresponds to the results shown in Ref. [9], 
in which we found quickly oscillating solutions with i? strongly deviating from Rgr- 

We consider a spherically symmetric cloud of matter with initially homogeneous energy density inside the limit 
radius r^- We choose the Schwarzschild-like isotropic coordinates in which the metric takes the form 

ds^ = Adt^ + B Sij dx^dx^ , (2.8) 

where A and B are functions of space and time. The corresponding expressions for the Christoffel symbols and 
Ricci tensor can be found in Appendix A. As usual the metric and the curvature tensors are expanded around their 
background values at first order in infinitesimal perturbations, i.e. 

A = Af) + 6A , 

B — B}j + 5B , ( 2 - 9 ) 

R = + 5R. 

Since Af, and R^ are quickly oscillating functions of time with possibly large amplitude (’’spikes”) which were found 
in our previous works and their time derivatives could be large, we shall keep terms of the second order in dt , such as 
d^A, dtAdtR, and so on. 

As it was done previously in our paper [11], we describe the evolution of perturbations using the equation for Gtt, 
the 9i9j-component of the equation for Gij, and the Euler and continuity equations which follow from the covariant 
conservation conditions = 0 and = 0. However, in the present work we do not conhne ourselves to a 

static background metric, but take into account its time variation which leads to several new effects. 


B. First order perturbations 


The tt-component of Eq. (2.4) is: 


1 


AH 


H2 3w2H 


AH- 


3Hi? d^Bd,R 


2A 


2B 


= Q, 


( 2 . 10 ) 


because, according to (A7) and (A14), Gtt = — AAH/H^ and Tu = qA. We take the background as weakly dependent 
on space coordinates but quickly oscillating with time, so that we can neglect the last term in the l.h.s. of Eq. (2.10). 
The corresponding equation for perturbations takes the form: 


where 5R is [see (A6)]: 


- MB -2Qb5B+^[ ASR - ^6BRb ) = Sg, 


SR = ASA - 3SB + 2ASB + -Af, SB. 


( 2 . 11 ) 


( 2 . 12 ) 


We assume that Hf, is small in comparison with Af, and Rb and put Af, = Hf, = 1 in denominators, see Eqs. (A17) 
and (A18). Analogously we find the 9^9,-component of the equation for : 


d,dj [SA + SB 


2SR\ _ 


3w2 J 


= 0 . 


(2.13) 


which remains unchanged with respect to Ref. [11]. The continuity equation, derived from = 0, is not changed 

either: 


Sg + gb djU^ + -gb SB = 0, 


( 2 . 14 ) 
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while the Euler equation acquires an additional term Af, gh Ujl2: 

9b j -\- djP +—£<{) djSA + 2 ^^ 9bUj = 0. ( 2 . 15 ) 

Such term obviously vanishes in static background, therefore it had been previously neglected. Notice that despite 
quick oscillations of Ah and Rb, the background energy density varies slowly with time, in agreement with the solution 
found in our paper [9]. 

Introducing = —Uj = —dja, P = and looking for the solution in the form ~ exp[—ik • x], we obtain the 
following system of equations for the five unknown functions of time, 5A, SB, SR, Sg, and cr: 


Sg= {e - 2gh)SB - SR-^SB, 
SR^^u}'^{SA + SB), 

< SR = -k^{SA + 2SB)-3SB + ^AbSB, 
Sg + gb a + -gb SB = Q, 
gb& - c^Sg + l-gb{Aba - SA) = 0 . 


(2.16a) 

(2.16b) 

(2.16c) 

(2.16d) 

(2.16e) 


At this stage we should comment on the application of the Fourier transformation to the derivation of equa¬ 
tions (2.16). The Fourier transform over space variables can be applied straightforwardly only if the coefficients in 
these equations are space-independent. It might be so in physically interesting cases, however in our example the 
metric function Ab depends explicitly upon space coordinates, in fact [see (A18)]: 

Ab = l + ^. (2.17) 

6 

Evidently, for sufficiently small r, such that Rbr'^/Q <C 1, so that we can completely neglect this term, the Fourier 
transformation can be safely applied. However, neglecting this term we exclude the effect of time dependent back¬ 
ground, which leads to new forms of instability. 

Nevertheless, we can take into account the effects of r-dependent terms and justify the applicability of the Fourier 
transformation both in the limits of large and small r. 

In the care of large r the Fourier transformation is applicable in the adiabatic case, when the background quantities 
change little over one wavelength of perturbations, A = 27r/fc. Indeed, performing the Fourier transform of the 
equations essentially consists in multiplying them by exp(ik • r) and integrating over cPr. If the coefficients of the 
linear equations for infinitesimal perturbations do not depend upon space coordinates we come as usual to an algebraic 
system of linear equations for the Fourier modes of the fluctuations. If, as is our case, some coefficients in the original 
differential equations depend upon r, then we can still transform our equations taking the integral exp(zk • r) not 
over the whole space but in some neighbourhood of a fixed value r = ro with radius Ar, chosen so that the space 
dependent coefficients are practically constant in such neighbourhood. To this end it is necessary that the range of 
integration contain many wave lengths, i.e. kAr ^ 1, and at the same time that Ar be sufficiently small, so that the 
r-dependent coefficients may be considered as approximately constant. Due to the largeness of kAr, such an integral 
would be close to the real Fourier transform with infinite integration limits. 

This essentially allows us to treat r as a constant parameter. Our final results will of course depend on such 
coordinate, in fact we will present results obtained in the two different limits fcr 3> 1 and /cr <C 1. The former 
condition is the more natural one, because the adiabatic approximation essentially implies kvm ^ 1 and most values 
of r are smaller but roughly of the same order of magnitude of the limit radius r^- The latter condition holds only 
near the centre of the cloud, that is for relatively small values of r; however, this may be a non negligible portion of 
the total system. The case of small r is considered below, after Eq. (2.24). 

Let us examine this condition more quantitatively. We demand that perturbations vary quickly compared to the 
background, that is 


drAb 



Ab 

5 

drAb 


(2.18) 


which, using (2.17), implies 


\Rb\r^ < 1 


kr 1 


(2.19) 
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Hence, choosing wavenumbers k ^ (see discussion above) and considering <C 1, as is the case for dilute 

enough systems, these conditions are nicely fulfilled and the Fourier expansion is reliable. The mass of the object 
under scrutiny is Mtot = 47rgbr^/3. We express Qb as 


Qb = 


STTgfc 

mh 



( 2 . 20 ) 


where Vg = 2Mtot/'m?pi is the gravitational (Schwarzschild) radius of the system. Correspondingly (2.19) becomes 


Qb irg ’ 


( 2 . 21 ) 


which sets an upper value for the amplitude of oscillations at which our approximations fail. Such constraint leaves 
a large portion of solutions as suitable for our analysis; for instance, for a solar mass cloud with mass density 10“^"' 
g/cTc? the ratio is Vmlfg ~ 10^^. 


C. Evolution equation 


Based on the comments above we derive from the system of five low order equations (2.16) the fourth order equation 
for the single function 5B\ 




Qb 






Ab 1 

' 2 4fc2 


1 + 


^2 J 


^ {^Rb + AbRb + 2RbC^k^^ 


Ab Ab , Ab 


fc2 


w" 1- 


3w2 J 

Qb 


Rb 


AbRb \ M, 


2fc2 


4^4 


SB 


SB 


clk^{k^ + w^) — 


4k^ 


3w2 


1 + — ‘^CgQbUJ 1 + —^ 


2k 


2 \ 1 


3w2 


SB = 0. 


( 2 . 22 ) 


In a wide parameter range, essentially whenever (2.19) holds, the following conditions are fulfilled: 

Al<^Ab, 

AbRb "C iib, 

AbAb <C Ab , 

AbRb "C iib, 

therefore several terms in Eq. (2.22) can be neglected. Using Eq. (2.17) we obtain: 


(2.23) 




4g& Rb 


1 + 


2k^ 


k'^r'^ 


SB + 


Rb 

(?gk^{k’^ + Cij^) — 




Qbi^ 


1 + 


Rbr^ 

12 

Ak?' 


Qb 




2*2 


- 2clQbi^J 1 + 


_ ^bcl ^ 2*2 \ 
3a;2 J 


SB + 


2*2 \ 


SB = Q. 


(2.24) 


Keeping both * and r in this equation looks eclectic. Nevertheless, it can be done in the adiabatic limit, as it is 
discussed above. Moreover, the limit of small kr makes sense as well. The applicability of this approximation can be 
justified as follows. One sees that Eq. (2.24) always contains the combination (1/*^ + ^2/6), where according to the 
arguments presented above, is understood as an adiabatic variable, so we can use both * and r. Adiabaticity is 
valid for kr ^ 1, so if we want to go to anti-adiabatic limit, we need to take the explicit Fourier transformation: 

(1/6) J d^rr'^e*'^'^SB{r) = -{l/6)dl[SB{k)] (2.25) 

and to compare it to SB{k)/k^. Evidently the latter dominates at small * if SB{k) is not singular at * —>■ 0. If SB{k) 
has a power law singularity, as Ijk'^, with v "^1, then numerically the 1/*^ term also dominates. We will arrive to a 
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similar conclusion working in the coordinate space, where we need to compare SB{r) r^/6 with 






For sufficiently small r the first term is again subdominant, so Eq. (2.24) is valid in the case of small kr. 
Next let us estimate the factor kr near the Jeans wavenumber k = kj = A/ph/(2Cg). We have 


f i, \2 3 rgr'^ 


Cs > i"g/r„ 


Even for a very cold gas at T = 10 K the speed of sound is about cl ^ T/rrip 


^ thus the above condition (2.27) 


can be fulfilled for appropriate densities. On the other hand, for extremely systems such as e.g. neutron stars the 
speed of sound is close to one (in units of the speed of light), and the ratio rgjrm is also much larger. However, in 
such systems the background pressure is of course non vanishing, while here we study non-relativistic systems with 
Ph « 0, so this is already outside of the limits of validity of our assumptions. 

We introduce the dimensionless time r = ujt, and define: 

6B = z , 

Bb — QbV : 


X = kr. 

Now we can rewrite (2.24) as the dimensionless equation: 


y'z'"+ 1 - 


+ 6(1 + c) + a 1 H—-—h 


+ ^(1 -c)]y' + 


be y' z' 


+ 6c(l + b) - 


— 2abc ( 1 + 


z = 0, 


where a prime denotes derivative with respect to r. 

We can further simplify the notation by introducing the following parameters: 


a = 1 + TT 


y- = b c(l + 6) - ^ f 1 + - 2ac f 1 + ^ 


(2.30a) 

(2.30b) 


fx = b c(l + 6) - - 1^1 + y j - 2ac (^1 + y j ■ (2.30c) 

We will make no assumption on the maximum value of k, hence on the maximum value of b which could be of order 
unity or even larger, but we will usually consider k > p/c^ and <C 1, corresponding to non-relativistic sound speed 
and wavenumbers at least of the order of the Jeans scale, for which GR predicts no exponential growth. This results 
in a <C 1, c <C 1, so in numerical calculations we keep only first order terms in these parameters. For the moment, we 
will put no restrictions on b and y. 

Let us rewrite (2.29) in the form convenient for qualitative and quantitative analysis: 

+y+(2a+ y"l z"+\(a+^) y'" - ^(1 + b)y] z'+ yiz = 0 . (2.31) 


Since the physically interesting quantity is the magnitude of density perturbations, we present 6q/Q b expressed through 
2 = 5B (e.g. in the limit of small x)' 


a{l + 2b/3) 


-2 +-z'y' + 


a(l -b 26/3) ■ 
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According to Eq. (2.20), when k is close to its Jeans value a = Qblk"^ ^ and the first term in the square brackets 
dominates if < 1 / 2 , which is surely the case for non-relativistic systems. 


III. MODIFIED JEANS INSTABILITY 


First of all, let us compare our equation (2.31) with the usual Jeans equation (remember the definition g = SirGg) 


Sg+(^c'lk'^ - Sg = 0. (3.1) 

As is well known, the sign of the term — gf ,/2 determines the stability of solutions, which are sound waves for 
> Qb/‘2' and unstable modes in the opposite regime. These unstable modes appear at the Jeans scale 


< ikrr s ||. (3.2) 

In modified gravity the condition of stability is determined by the sign of g, (2.30c). In the limit of small amplitude of 
curvature oscillations or very low frequency, we can neglect y(r) in Eq. (2.31) so that it is reduced to a simple equation 
with constant coefficients which is solved by the substitution z = exp( 7 T). The eigenvalue 7 is thus determined by 
the algebraic equation: 

7 ^ T 11^7^ “b M ~ 0 , (3-3) 

where 11^ is given by Eq. (2.30b). The eigenvalues 7 ^ solving this are: 

7 =-—±\l—-g. (3.4) 


If ^ < 0, then one of the roots 7 ^ > 0, so one of the eigenvalues is positive. It corresponds to the usual exponential 
Jeans instability, though the values of the Jeans wave vector in modified gravity and in GR are different. The 
magnitude of the Jeans wave number is found from the condition g = 0, which in the case of small sound speed 
(Cg <C 1 ) yields: 


2 c(l + b) 
~ 1 + 45/3 ■ 


(3.5) 


This is an equation quadratic with respect to the Jeans wave number [see the definitions (2.28)]. We present an 
explicit solution for large ^ k"^, g: 


r /r,GR',2-| 

(fc“G)2 ^ (^GR)2 1 + i^ , 


(3.6) 


which recovers the GR result in the limit w —>■ 00 . For large w, the physical eigenvalues corresponding to (3.4) are 


t ^2 22 

^mg — w 7 _ 


Qb 


— cik'^ + 


— + 


0b 

y 

gbbo" 


Qb 

2k^ 


— cz 


k^ 


Qb 




(3.7) 


The former value gives the usual Jeans growth, up to a small correction (the term in brackets), whereas the latter 
corresponds to scalaron oscillations with frequency approximately equal to w, as we expected since the theory possesses 
an additional massive scalar mode. 

Equation (3.6) shows that in modified gravity the Jeans wavenumber is larger than in GR. This corresponds to a 
reduced minimum length scale associated to structure formation. The correction is typically small, but for models 
in which fcj^/w is non negligible it could lead to significant corrections. The phenomenology of this result and its 
implications for constraining models will be dealt with elsewhere. 

If g is positive, but g < 11^/4, both possible values of 7 ^ are real and negative, so 7 is purely imaginary which 
corresponds to acoustic oscillations. Thus these two cases of negative and positive (but not too large, see below) g are 
in a one-to-one correspondence to the usual Jeans analysis. The results of this paper generalize those of [II], where 
we analyzed gravitational instability in modified gravity over a time-independent background. 



If /i is large, i.e. /i > n‘^/4, there would exist a new type of unstable oscillating solution with exponentially rising 
amplitude. Indeed, 7 ^ becomes a complex number and two of the four solutions of Eq. (3.3) for 7 have positive real 
parts. However, such solution is absent in the model studied here; using (2.30), we find that the condition p. > 
becomes: 


(4a + 3c)^ 


+ l-2c 


b^- 


a(5 — 4a — 21c) 


— 2 -t- 2c 




(3.8) 


This equation has a real positive solution for b if and only if the coefficient [a(5 — 4a — 21c)/3 — 2 + 2c] is positive. 
However, this is not fulfilled for any positive a and c. On the other hand, it is not excluded that some other modified 
gravity models may possess such a kind of instability. 


IV. PARAMETRIC RESONANCE 


In this section we show that even for /x > 0, corresponding classically to the stable (sound-wave) regime, Eq. (2.31) 
may admit unstable solutions, through a mechanism analogous to parametric resonance. 

The usual textbook example of parametric resonance is the Mathieu equation (for details, see e.g. [16]): 

/(r) +nl[l + h cos(Hi -h e)T] /(r) = 0 , (4.1) 

where it is assumed that h <C 1. If Hi/Hq = 2 , the solution grows exponentially, behaving approximately as: 


with^ 


Parametric resonance is excited for 


/(r) - sm(Hor -h (/?) exp( 7 T), 


7 = 


h 

~T~ 


kl< 


h 

~Y~ 


(4.2) 


(4.3) 


(4.4) 


There are other resonance modes, which are generated at the frequencies Hi/Hq = 2ln, where n is an integer, but 
they are usually noticeably weaker. 

In modified gravity, we have derived the 4th order equation (2.31) governing the evolution of scalar perturbations 
over an oscillating background, which enters through the function j/(t). An effect analogous to parametric resonance 
appears in this situation as well. A similarity of parametric resonance in the fourth-order equation and in the classical 
Mathieu equation can be seen both in the limit of high and low frequencies. In the former case one can neglect lower 
order derivatives in Eq. (2.31), keeping only z"", z'", and z", effectively reducing the equation to a second order one 
in the variable z". In the case of low frequencies, on the other hand, z"" and z'" can be neglected and we come 
again to the Mathieu equation. In what follows we study the complete fourth order equation both analytically and 
numerically and observe resonant amplification. 

We assume that y(r) is a periodic function describing curvature oscillations and consider two possible forms of 
them: 


1. purely harmonic ones (see Sec. IV A): 


Uharmit) = yeq{t) -|- j/Q cos(a;it -|- 9 ), (4.5) 

where yeq is the equilibrium point of the potential around which the curvature oscillates, i.e. the external energy 
density in units of the initial density yeq = git)/go. The harmonic solution ytarm was found in our paper [9] 
when the amplitude of curvature oscillations is small. In this case uji = uj and yo < yeq- The time variation of 
yeq and yo is much slower than the oscillations, namely yeq/yeq ‘C uji and similarly for yo. 


^ The parameter 7 introduced here is evidently different from the growth rate 7 characterising the Jeans instability, considered in the 
previous section. 
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2. spiky solutions found in Ref. [9] (see Sec. IV B), which, up to a slowly changing term, we approximate as 


ysp{t) 


yo (P 

(P + sin^(a;2t + 9) ’ 


(4.6) 


where d <C 1, so that we have spiky solutions (narrow peaks with a large separation between them). 

According to [9] we have W 2 = w/2, where uj is given by Eqs. (2.5) and (2.6). The Fourier transform of ysp{t) 
contains modes with much higher frequencies than u) 2 - The dominant mode, i.e. the mode with the largest 
amplitude, is excited at 2 w 2 . 

When studying the spiky solutions, and in general any non-harmonic behaviour of the background curvature, we 
will have several modes with frequency at least equal to w and possibly much larger. The overall effect can be 
seen as the combination of the individual contributions from the different modes, so it is useful to keep wi or W 2 
as a free parameter. If resonance is excited for some value of a;i ,2 and some non-harmonic solution for R contains 
such frequency, then we can expect resonant behaviour in the full non-harmonic case too, though perhaps slightly 
suppressed. 


A. Harmonic Oscillations 

Let us start from harmonic oscillations, Eq. (4.5). We assume that curvature oscillates with a fixed frequency. 
Indeed, this is a rather general scenario in F{R) theories, at least when the conditions (2.2) and <C I hold.^ In 
terms of dimensionless time, r = wt, yharmir) takes the form: 

2 /harm( t) = yeq{T) -f 7/0 COs{^iT ^ 9) , (4.7) 

where = wi/w, t/o is the amplitude of oscillations and d is a constant phase. This form of y physically corresponds 
to an oscillating scalar curvature [see (2.28)]: 

Rb{t) = -gb[yeq{t) + 2/0 cos(wit -b 9)] = RoRiyeq + 2/0 cos(a;it -b 9)] . (4.8) 

In order for this to be physically sensible we need 2/0 < 2/eg 7 so that R does not change sign, but note that T/eg rises 
with time because the background energy density is assumed to be constantly increasing, see e.g. the second paper 
in Ref. [8] and Ref. [9]. Therefore, it makes sense to explore values 2/0 > 1- 

Eq. (2.31) is reduced to the Mathieu equation (4.1), if the terms containing odd derivatives in Eq. (2.31) can be 
neglected and /i <C 1. We can see that this is the case when a, a, b and c are all much smaller than unity. Although 
this seems rather restrictive, it is actually a natural possibility as it corresponds to A: <C w, <C 1 and k'^ ^ Qb- In 
this case the substitution z" = x yields: 

x" -b -b 2a cos(nir)] a; = 0 , (4.9) 

If Oi/n = 2, then the parametric resonance would be excited. Comparing this to Eq. (4.1), we see that the parameter 
h is expressed through a as 

h = 2cAI“2 . (4.10) 


1. Fundamental harmonic 


Let us consider = 1, as is the case for harmonic oscillations with small amplitude, in which the curvature 
oscillates with a frequency given precisely by (2.5). Then the parametric resonance condition requires 



1 - 


b- 


Aab 


1 - 


Qb 


^- 2~ 


^Qb 
3^2 ’ 


(4.11) 


® If the frequency changes with time, the resonance effect may possibly survive but with a changing frequency of the signal as well. 
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where we express ft in terms of physical quantities, according to Eqs. (2.28), (2.30a), and (2.30b) assuming c <C 1. 

Due to the similarity to the standard Mathieu equation, we are tempted to study how our full equation (2.31) 
behaves when (4.11) is fulfilled. Solving for k gives 

= ^( 16p& - 9^2 ± + 256^ ) ~ i (4.12) 

3 

The former solution is of course unphysical because it leads to an imaginary^ k. The latter solution 

kL ^ ^ (4.13) 

is physically sensible but lies outside the Jeans radius, where we have unstable modes even in GR; still, provided that 
k-cesfm > 1 so that we are still inside the cloud considered, this would result in an enhancement of structure formation 
starting at the specific scale (4.13). The implications of this will be considered elsewhere. 


2. Higher harmonics 


As mentioned earlier, when curvature does not behave precisely as (4.7) due to anharmonic features, we will have 
several modes contributing to the overall effect, and these modes will in general contain higher frequencies. Therefore, 
let us now keep Di as a free parameter and discuss parametric resonance in this case. 

As in the standard case, we assume that Eq. (4.7) is fulfilled and insert the tentative solution 

z{t) ~ (A cos Dr + R sin Dr) exp(yr) (4-14) 

into Eq. (2.31). The transformed equation would contain higher harmonics, which we neglect (as is done in the 
standard case). In order to do so, we exploit the following exact relations: 


2 cos X cos 2x = cos x + cos 3x, 

2 cos X sin 2x = sin x + sin 3a;, 

2 sin X cos 2x = — sin x + sin 3a;, 
2 sin X sin 2a; = cos x — cos 3x , 


(4.15) 


and neglect cos3x, sin3x. We also take the limit 7 <C D, as expected and confirmed by results (see below). We 
will investigate values Di/D ~ 2, which as we have seen above is equivalent to the classical condition for parametric 
resonance. 

We are left with an equation containing terms proportional to cos Dr and sin Dr. We demand that both coefficient 
vanish simultaneously, which yields: 


(Pi+P2)A + 7(Qi+Q2)S = 0 
7(Qi-Q2)A + (Pi-P2)P = 0, 


where 

Pi = 24/u., Qi = -j/oD[8D^(3a + ay^) + (a - 6 q;)x^] , 

P 2 = yofi^Kax^ - 6a(x^ + 4D^))] Q 2 = -48D^ . 

This system admits non trivial solutions when the determinant of the associated matrix vanishes, namely when 


(4.17) 


pJ _ pz 

72 ^ ^2 


Qi Q 2 


(4.18) 


4 


In inhomogeneous systems of finite size k may be imaginary, as is known for propagation in waveguides. 
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The definitions (2.28, 2.30) in the limit of large oj yield the following result for the physical growth rate in the resonant 
case: 


1^2 2 2 
^res OJ 


l/§gg(6 + X^) 

576fc4 


[(6 + x^)w^ + (14 + X^)fc^; 


2-^2 2 

yoQb^ 


2~2 2 4 

ypgfc^ X 
576*4 


(x < 1), 
(x » 1) ■ 


(4.19) 


The ordinary GR eigenvalue (3.1), which at the lowest perturbation order coincides with the modified gravity (non¬ 
resonant) solution (3.7), is 


-p 2 ■p 2 

^ GR — MG 




(4.20) 


The resonant behaviour thus dominates for 


Tgr 

T res 


< 1 


( Sk'^iQb - ‘2-clP) 




2/0 >< 


288(p 6 - 2c2fc2) 


^^2x4 


(x< 1), 
(x » 1) • 


(4.21) 


When studying the range of scales where we would have structure growth in GR (* < fcj), the factor [gb — 
is positive. In the opposite region, every resonant behaviour is of course the dominant one since GR simply predicts 
sound waves. 

Equation (4.21) indicates that there exists a lower limit for j/q to activate parametric resonance, as it happens in 
the usual case. Nevertheless, there should be no problem in finding a sufficiently large portion of the viable parameter 
space in which this condition is satisfied. 


3. Numerical Results 


a. Small radii. 

The condition x = ^ 1 corresponds to the innermost portion of the cloud, at distances from the centre shorter 

than Although this is likely a small region, it may still be non-negligible, especially when studying structure 

growth at relatively large scales (that is for relatively small k). 

We solved Eq. (2.31) numerically for different values of t/o and fli. The parametric resonance excitation is observed 
at the expected frequency rii/n = 2, see Fig. 1. We present the numerical solutions for z{t) and compare them with 
our analytical estimate for the amplitude of z, namely 

z--exp(rrest) , (4.22) 


with Ei-es given by (4.19). 

In this paragraphs the parameters of the medium and the wave number of the fluctuations were taken to be 

a = 6 = x = 0.01, c = 0.02, (4.23) 

and varying j/o and fli. These values were merely chosen to produce the figures showing the resonant behaviour, which 
however appears for a very wide range of the parameters. On the other hand, parametric resonance is rather sensitive 
to variations of j/o and of course of Oi. In fact, the exponential growth becomes slower when we move from Oi/fl = 2 
towards fli/n = 2.015 as one can see comparing the top and bottom panels in Fig. 1. Such sharp frequency dependence 
clearly demonstrates the resonance behaviour, with the resonance half-width roughly of the order of ~ 0.1 fl. Such 
a behaviour helps to distinguish between the parametric resonance instability and antifriction instability considered 
below (Sec. V), for which the frequency dependence is very weak. We have not observed resonance at fli = 17, 
expected in the Mathieu equation. 

In both cases considered in Fig. 1, the agreement of numerical calculations with the analytic estimate (4.19) is 
remarkable, and as good as iFMc/rnum — 1| ~ 5 x 10“^ for fli/n = 2 and decreasing to about 24% for = 2.015, 
which is close to the the resonance threshold. 
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FIG. 1. Solution for the parameters (4.23). Left panels'. Oscillations of « as a function of dimensionless time r. Right panels: 
\z{t)\ (note the logarithmic scale) and, in orange, the envelope (amplitude) of the oscillations, namely 4.22) with growth rate 
given by (4.19). 

Top: Parametric resonance excitation of z{t) for harmonic curvature oscillations with yo = 5 and fli/fl = 2 and comparison 
with the predicted exponential rise. The relative difference between the exact and estimated growth rates is about 5 x 10“®. 
Bottom: Results for fli/fl = 2.015, other parameters unvaried. Clearly the resonant behaviour is much weaker, but the agree¬ 
ment between the exact numerical value of P and the analytical estimate (orange) is still satisfactory: IPres/rnum — 1| 0.24. 

Note that although solutions are only shown for relatively large times, the evolution of the system starts at r = 0. 


b. Large radii. 

The case of large x or r ^ corresponds to a larger volume of the collapsing cloud. According to our assumptions, 
perturbations normally vary on scales much shorter than the total size of the object hence krm ^1. If we take 
e.g. 0.5 Cm < r < r^, this region would occupy about 0.9 of the cloud volume, and if 0.1 < r < r™, then the 

corresponding volume is 0.999 of the total volume of the cloud. 

We solved Eq. (2.31) numerically, and present our results for the following parameters: 

a = 6 = c = 0.01, x = 10, 2/0 = 1, (4.24) 


and varying fli. 

In Fig. 2 we show results for fli/Sl = 2 and fli/Lt = 2.075. As we move away from the resonant value of the 
frequency, the agreement with our analytical estimate (4.19) decreases but remains fairly satisfactory, the relative 
discrepancy being much smaller than an order of magnitude. The width of the resonance region for is roughly 
i5ni/r2i ~ 0.05. However, this depends strongly on the value of yo: for instance, for 2/0=4 the width is huge, about 
i5ni/r2i ~ 1.2. The resonance leads to a much faster growth than in the case of small x- 

Of course huge values of the amplitude of oscillations, which can be easily reached numerically, do not make much 
sense because the first order approximation used in deriving (2.16) are valid only if perturbations are much smaller 
than unity. With an initial amplitude of fluctuations of the order of 10“"'^ — 10“® (as from CMB data) the results can 
be trusted up to amplifications not exceeding roughly 5 orders of magnitude. 


B. Spike-like oscillations 

As it was found in our previous works [9], harmonic oscillations do not cover the entirety of possible solutions for 
the background R{t). In fact, we have shown the possibility of narrow “spikes” of large amplitude (|i?| ^ |i?GR|), see 
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FIG. 2. Solution for the parameters (4.24). Left and right panels depict the same quantities as in Fig. 1. Top: Solution for 
1/0 = 1 and fli/fl = 2 and comparison with the predicted behaviour, in orange. The relative difference between the exact and 
estimated growth rates is ~ 3 x 10“^. Bottom-. Results for = 2.075, other parameters unvaried. The resonant behaviour 

is weaker but still rather strong, but the agreement between the exact numerical value of F and the analytical estimate (orange) 
has decreased: IFres/Fnum — 1| ~ 0.32. Still, our analytical estimates are well within one order of magnitude from the numerical 
results. 


Eq. (4.6). These spikes have time separation equal to 27r[a;(i? = i?GR)]~^ and much smaller width; for instance, in 
the case of the model (1.2), the width of the spikes can be as small as roughly < (10^ GeV)“^. 

We solve numerically Eq. (2.31) using ysp as function of the dimensionless time: 


Vspir) 


yo<P _ 

d? + sin^(n2T + 6) ' 


(4.25) 


with 172 — 1/2 in agreement with the solution of Ref. [9]. This function and its Fourier transform are presented in 
Fig. 3. Note that the amplitudes of even harmonics are much larger than the amplitudes of odd ones. Moreover, the 
amplitudes of odd harmonics can be negative despite the fact that the curvature ?/(t) remains positive, as it should. 


.'<'■) yia,) 




FIG. 3. Left panel: y{r) dehued by Eq. (4.25) with yo = 1, 172 = 0.5, and d = 0.1. Right panel: Fourier transform of i/(t). 

Fig. 4 (left panel) clearly shows parametric resonance at 172 = — 0.5013 which corresponds to the second 

resonant mode, see Eq. (2.30b), with the chosen values of parameters: a = b = 0.01 and c = 0.02. For these values fi 
is positive and equal to 1.47 • 10“^. In the right panel z{t) is presented for a lower frequency 172/17 = 0.496. We see 
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that the resonance is still excited but much weaker. The shape of the resonance curve is symmetric with respect to 
the position of the resonance frequency. The sharp decrease of the oscillation amplitude shows that the resonance is 
quite narrow, having width Sfl/fl <C 1. 



FIG. 4. Parametric resonance excitation of z{t) for spike-like curvature oscillations with yo = 30 and = 1/2 (left panel) 

and 0 . 2 /^ = 0.496 (right panel). 

The main mode of parametric resonance should be at ^ 2 /^ = 1; it is presented in the left panel of Fig. 5. However, 
it appears to be sub-dominant with respect to the mode at f22/H = 1/2 because of the suppression of the odd Fourier 
amplitudes of y(T), see Fig. 3. 

In Fig. 6 the evolution of z(t) is depicted for H 2 /H = 1/3 (left) and H 2 /H = 1/4 (right). These higher modes are 
weaker as expected, and as mentioned above the Fourier amplitudes are weaker for odd modes. We see that the 4th 
order equation demonstrates parametric resonance effects quite similar to the classical case (4.1), though the impact 
of the odd derivatives in Eq. (2.31) may be significant and lead to quantitative modifications of the results. 



FIG. 5. Parametric resonance excitation of z(t) for spike-like curvature oscillations with yo = 30 and ft 2 /ft = 1 (left panel) 
and fl 2 /fl = 0.992 (right panel). 




FIG. 6. Parametric resonance excitation of z(t) for spike-like curvature oscillations with yo — 30 and 0,2/0. = 1/3 (left panel) 
and O 2 /O — 1/4 (right panel). 
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V. ANTIFRICTION AMPLIFICATION 

Antifriction amplification of perturbations is induced by the change in the signs of the coefficients in front of the odd 
derivatives in the 4th order equation. If they become negative (and they do), then instead of damping the oscillations 
they would lead to their enhancement, as can be easily seen solving the oscillator equation; 

z + Fi + = 0 (5-1) 

with r < 0. We see below that for sufficiently large a (2.30a) and/or y (2.28) the fourth order equation (2.29) has 
indeed strongly rising solutions. 


A. Harmonic curvature oscillations 

1. Large amplitude: analytic solution 

In section III we considered the case of small amplitude of the oscillating background and saw that equation (2.31) 
can be solved analytically. It is interesting that this equation can be also solved analytically in the opposite limit of 
large amplitude of oscillations, but with small c and y. In this case Eq. (2.31) becomes 

z'"' + ay'z"' + 2az"y" + az'y'" = z"" + a{z'y')" = 0 . (5.2) 

The equation is easily integrated: 

z" + az'y = Cl + C 2 T , (5.3) 


leading to the solution: 


z'= Co + Cl [ dn + C2 f dririe“^(^^) , 

Jo Jo 


or equivalently 


z = zo 


f dTi ^ dr 2 {Ci + C 2 T 2 ) e 


ay{T2) 


(5.4) 


(5.5) 


We take for definiteness ?/(t) = yQCOs{iliT). From the expression (5.4) it is clear that the derivative z' is small 
when ay < 0, while z' is positive and large for ay > 0, so the function z remains constant or grows in the first and 
in the second case respectively. 

This behaviour is shown in Figs. 7, 8. The amplihcation of z(t) takes place independently of the frequency of the 
background curvature oscillations wi. This can be explained by the change of the sign of the coefficients in front of the 
odd derivatives in Eqs. (2.31) and (5.2) for large yo- When these coefficients are positive they act as a friction force but 
when they are negative they act as antifriction. Such effect is showed numerically below for the full equation (2.31). 




FIG. 7. The function z(r) {right panel) and its derivative z'{t) {left panel), as given by Eq. (5.4) for Co = 1, Ci = C 2 = 0, 
j/(r) = yo cos(r), and ayo = 2. 
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FIG. 8. The function z{t) {right panel) and its derivative z'{t) {left panel), as given by Eq. (5.4) for C 2 = 1, Co = Gi = 0, 
?/(''") = yo cos(r), and ago = 2. 


2. Numerical solutions 


a. Small radii 

Here we solve eq. (2.31) numerically for the ad-hoc chosen parameters a = 0.01, b = 0.01, c = 0.02, x = 0.01, for 
different values of yo and Hi. The antifriction amplification is observed at the frequencies away from the resonance 
values, if yo exceeds a threshold value, yth- Tbe farther away the frequency is from the resonance, the larger the 
threshold. For example for uji = 3.2 the threshold value is yth = 169, while for wi = 4.4 the threshold is yth = 267. 
The evolution of z{t) in these two cases is depicted in Fig. 9. 


Z 




FIG. 9. Antifriction effect in the evolntion of z{t) for harmonic curvature oscillations with yo = 169, Hi/H = 3.2 {left panel) 
and yo = 267, H 2 /H = 4.4 {right panel). Small radii: x ~ 0.01. 


b. Large radii. 

A large x leads to an antifriction instability more efficiently, because the term proportional to enters only into the 
coefficients in front of the first derivative, z', and there is no destructive interference of the first and third derivative 
terms, which in the general case may have opposite signs and act in opposite directions. 

In Fig. 10 we present an example of antifiriction behaviour in systems with harmonic curvature oscillations for the 
same values of parameters a = 0.01, b = 0.01, c = 0.02 and frequencies as in the case of small radii, except for the 
value of X- We take x = 10 and again observe the threshold behaviour, but for the frequency Hi = 3.2 the threshold 
amplitude is yth = 10, whereas for Hi = 4.4 the threshold amplitude yth = H- 


B. Spike-like curvature oscillations: numerical solutions 

1. Small radii 

We present here the numerical solution of Eq. (2.31) for spike-like excitations (4.6) with different frequencies H 2 
and amplitudes yo. We see that for all frequencies, not only the resonant ones, the solutions are unstable, provided 
that yo is larger than some threshold value which depends on the frequency. We attribute this phenomenon to a 
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FIG. 10. Antifriction effect in the evolution of z(t) for harmonic curvature oscillations with i/o = 10, fli/fl = 3.2 {left panel) 
and yo ~ 11, = 4.4 {right panel). Large radii: x ~ 10- 


change of the sign of coefficients in front of the odd derivatives in Eq. (2.31). In a sense this result is similar to that 
described in the previous subsection for harmonic j/(t), though quantitatively different, both in magnitude and in the 
shape of the signal. The latter may be explained by the different form of j/(t), by different initial conditions, or by 
other effects induced by a non-zero ^ neglected in Eq. (5.2). 

In Fig. 11 we present the evolution of z{t) for out-of-resonance frequencies Q. 2 /^ = 0.6, 0.7. The farther away 
the frequency is from the resonant one, approximately equal to 0.5, the larger the threshold value of j/o necessary for 
generating an unstable solution. The magnitudes of ?/o taken in these figures are quite close to the threshold values, 
2/0 = 400 for 0 , 2 /^ = 0.6 and yo = 702 for 0 , 2 /^ = 0.7. Note that the threshold value of yo at the resonance with 
122/^4 « 0.5 is about 15. This threshold effect is also present in the standard parametric resonance phenomenon, if 
friction is non negligible [16]. 



FIG. 11. Antifriction effect in the evolution of z{t) for spike-like curvature oscillations with yo = 400, Q. 2 /O. = 0.6 {left panel) 
and with yo = 702, Q. 2 /O. = 0.7 {right panel). 

In Fig. 12 the evolution of y(r) is depicted for higher frequencies il 2 /^ = 0.8 {left panel) and D. 2 /^ = 0.9 {right 
panel) with yo = 1280 and yo = 372 respectively. The higher frequency ^ 2 /^ = 0.9 is closer to the resonant one 
122 ~ so the threshold value of yo is smaller. 


2. Large radii 

The results of numerical calculations are presented in Fig. 13. The instability is very clearly seen, but the character 
of the two kinds of instability for different values of the parameters, depicted in the left and right panels are very 
much different. In the first case we observe oscillations with quickly increasing amplitude, while in the second case 
there is an explosive, practically monotonic rise. 
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FIG. 12. Antifriction effect in the evolution of z{t) for spike-like curvature oscillations with yo = 1280, 0,2/0. = 0.8 {left panel) 
and yo = 372, O 2 /O = 0.9 {right panel). 



FIG. 13. Antifriction effect in the evolution of z{t) for spike-like curvature oscillations with yo = 14, fli/fl = 0.4 {left panel) 
and yo = 9, O 2 /O — 0.6 {right panel). Large radii: x = 10- 


VI. DISCUSSION AND CONCLUSIONS 


A general feature of F{R) modified gravity are high frequency oscillations of curvature and metric in contracting 
matter systems such as a gas cloud of interstellar or intergalactic matter in the process of star or galaxy formation. 
A similar effect might take place in the process of stellar collapse, e.g. before a supernova explosion and subsequent 
formation of neutron star or black hole. As we have shown in this work, the evolution of metric and density per¬ 
turbations in an oscillating gravitational background possesses new and very interesting features in addition to the 
usual exponential Jeans instability, which itself occurs at a different scale than in GR. The oscillating background 
metric and curvature of the system induce an effect analogous to the standard parametric resonance, resulting in 
an amplification of the fluctuations. There is a quantitative difference with respect to the standard case because 
the usual parametric resonance is described by a second order differential equation, while in modified gravity the 
equation (2.22), which governs the evolution of fluctuations, is fourth order. We study here a rather general equation 
without restrictions to any specific J^(i?)-theory. Still, we used as a guiding example the model described by the 
action (1.2) but did not confine ourselves to specific values of the parameters demanded by this model. We used a 
similar approach concerning the form of the oscillating metric and curvature background, considering the two cases 
of arbitrary harmonic oscillations and spike-like solutions of the type found in our previous work [9]. In both cases 
the evolution of perturbations is described by a fourth order equation, in which odd derivatives may play a crucial 
role. They describe the damping due to friction if their coefficients are positive (this is the usual case in the standard 
parametric resonance theory). However, such coefficients may periodically change sign and this leads to a periodic 
antifriction force and to a consequent amplification of perturbations. The latter effect occurs at rather high amplitude 
of curvature oscillations, higher than the background GR value. Systems where curvature oscillations with large 
amplitude are induced have been found in our work [9]. 

The fourth order equations (2.29) or (2.31), which govern the evolution of the density perturbations, demonstrate 
a very rich pattern of different types of instabilities. There is a close analogue of parametric resonance, which is 
easy to describe theoretically, almost in the same way as the usual parametric resonance. However, in addition to 
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FIG. 14. Evolution of &qIq in parametric resonance region induced by the spike-like solution with yo = 30, Q. 2 /Q. = 0.5. The 
initial value of metric perturbation is taken as 5B{Q) = 2(0) = 10“®. 


this we found a new kind of instability induced by negative signs of the coefficients in front of the odd derivatives 
in the equation - we dub it the antifriction instability. Interestingly this instability generates a completely different 
behaviour of the rising perturbations, depending upon the values of the system parameters. There can be an oscillating 
behaviour with quickly rising amplitude, or a quasi-explosive one with amplitude tending monotonically to infinity. 
Of course these results are reliable only if perturbations are sufficiently small, not larger than unity when written in 
dimensionless form. 

We have studied (2.29) analytically and numerically, solving for the perturbation in the space-space metric term. 
Based on these solutions, one can calculate the evolution of the relative density contrast using (2.32). Barring 
accidental cancellations, exponentially growing solutions for z will lead to an equivalent behaviour for the density 
perturbation 5q, as shows in Fig. 14 in the case of parametric resonance induced by the spiky solution with moderately 
large amplitude yo = 30. 

We can see that in a relatively short time, about 10^0;“^, the density perturbation would reach unity if the initial 
metric perturbations are about 10“^. The effective time can be much shorter than the gravitational time typical for 
the Jeans-type rise of density perturbations, since w is usually large. Though the density contrast is an oscillating 
function of time, its impact on structure formation may be non-negligible and should lead to constraints on the 
parameters of the underlying F"(i?)-theory. This will be studied elsewhere. 

Another potentially interesting effect of the oscillating rise of curvature perturbations is the generation of curvature 
propagation around contracting bodies, i.e. a scalar mode for gravitational waves having mass equal to w, which may 
vary from practically zero up to the scalaron mass m. They may lead to a considerable loss of energy of the collapsing 
objects. 
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Appendix A: Metric and Curvature 


As we did previously, we consider a spherically symmetric cloud of matter with initially constant energy density 
inside the limit radius r = Vm- We choose Schwarzschild-like isotropic coordinates in which the metric takes the form: 

ds^ = Adt^ — B Sij dx^dx^ , (Al) 


where the functions A and B may depend upon r and t. The corresponding Christoffel symbols are: 

t A , djA , pf^dkA , SjkB 

Itt 2A ’ 2A ’ 2B ’ 2A ’ 

6^ B 1 

^ + 8piB - %J'="5„R). 


(A2) 
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For the Ricci tensor, including terms quadratic in F’s, we obtain: 

A A 3B 3^2 3AB d^AdjB d^AdjA 

4^ 4AB 4R2 4AB ’ 

d,B , Bd,B , Bd,A 
B B^ 2AB ’ 

^ ( B AB B"^ AB d'^AdkB ^ d’^BdkB\ 

% - V + 4AB “442 4AB ^ 4 R 2 j 

d,djA d^djB , ^ Sd^BdjB ^ d^AdjB + d^Ad.B 

2I 2^ 4^2 452 + 4AB ■ 

Here and in what follows the upper space indices are raised with the Kronecker delta, A = S^^dkA. 
The corresponding curvature scalar is: 

„ AH 3H 2AB iAB d^Ad,A Bd.B d^AdiB 


^ AB “ AB B2 2A^B ~ 2A^B 2B3 
Let us now present the expressions for the Einstein tensor Gfj^, = R^i, — \/2g^i,R-. 

_ AAB 3B2 iAd^BdjB 

tt - ^ ^ + ^^3 ’ 

^tj ^tj 5 


2 HB 2 ■ 


^ ^ (AA ^ AB B ^ B^ ^ AB d'^AdkA d'^BdkB] 

~ W “ H 4AB ^ 442 2B2 j 

didjA didjB diAdjA SdtBdjB diAdjB + djAdiB 

2 I W~ 4H2 4B2 4AB ’ ' 

The energy-momentum tensor is taken in the perfect fluid form without dissipative corrections: 

T^u = {q + P)Ufj,Ui/— Pgfiy ^ [Jr 

where g and P are respectively the energy density and pressure of the fluid and the four-velocity is: 

C/^ = — and C/^= 5 ^„C/“. (; 

We assume that the three-velocity = dx^/dt is small and thus neglect terms quadratic in v. Correspondingly, 

U- = _^ u 

^ \fA^\ — {BlA)vjvi VA 

According to our definition Vj = . From the condition 


1 = g'^'^U^U, = jU^ - ^6’^^UkU, « jU^ 


we find Ut ~ VA. Now we can write: 

Ttt = {g + P)Uf -PA^gA, (A14) 

Tjt = ig + P)UtUj --{g + P)Bvj , (A15) 

Pj = (p + P)U^UJ - Pgrj « PB5^J. (A16) 

Background metric for spherically symmetric distribution of matter (an analog of the Schwarzschild solution for 
modified gravity) has been found in several works. We use here the form for the internal solution obtained in our 
paper [13] (references to other papers can be found there): 




A,{r,t) = l + ^^+Af^'^\rA), 
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where 


M{r,t)= [ (frToo{r,t)=4Trl' dr Too{r,t), 
Jo Jo 




3r„ 


'^Qm ^^2 

2r3 2r„, 3m 


,2 

'-Pi 


f 2 2\2 

(»'m - ry , 


and Tg = 2M/m^pi with M being the total mass of the object under scrutiny. 


(A19) 

(A20) 
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